function [elements,processInfo] = CalculateMaterialResponse...
    (elements,model,processInfo)
ne = model.countElements();
for i=1:ne
    e = model.getElement(i-1);
    f = e.getFace();
    ue = f.getUHatLocal();
    w = e.getGaussWeight(0);
    xi = e.getGaussPoint(0);
    B = zeros(3,6);
    B = e.computeB(f,B,xi);
    C = elements(i).constitutiveMatrix;
    if elements(i).isCrack==1
        jac = f.jacXAt(xi);
        detJac = det(jac);
        %update Xi
        due = ue - elements(i).ue_old;
        dxi = elements(i).Kxixi\(elements(i).renh-...
            elements(i).Kxid*due);
        elements(i).Xi = elements(i).Xi+dxi;
        %compute strain
        epsilon = B*ue + elements(i).Gc*elements(i).Xi;
        %compute stress
        elements(i).stress = C*epsilon;
        %compute TGamma and update Kxixi
        lGamma = elements(i).lGamma;
        J = elements(i).J;
        Xi = elements(i).Xi;
        u_mu_n = Xi(1);
        t = elements(i).ft+elements(i).H*u_mu_n;
        TGamma = [max(t,0) 0]';
        if t < 0
            %fully opening crack
            elements(i).Kxixi = elements(i).Ge'*...
                C*elements(i).Gc*(w*detJac);
        else
            %partial opening crack
            Ctilde = [elements(i).H 0;0 0]*J';
            elements(i).Kxixi = elements(i).Ge'*...
                C*elements(i).Gc*(w*detJac) + ...
                lGamma*J'*Ctilde*J;
        end
        %compute renh
        renh = -lGamma*J'*TGamma - elements(i).Ge'*...
            elements(i).stress*(w*detJac);
        elements(i).renh = renh;
    else
        %compute stress
        elements(i).stress = C*B*ue;
    end
    %save ue
    elements(i).ue_old = ue;
end
end
